library(tidyverse)

q <- seq(.01,.99, .01)

p00 <- .9650
p11 <- .0016
p10 <- .0202
p01 <- .0132

id_min <- id_max <- rep(NA, length(q))

for (i in 1:length(q)) {
  id_min[i] <- max(-1*(p10 + p11)/(1-q[i]), 
                   -1*(p01/q[i])-(p10/(1-q[i])),
                   -1*(p00 + p01)/q[i])
  id_max[i] <- min((p00/(1-q[i]))-(p01/q[i]),
                   (p11/q[i])-(p10/(1-q[i])))
}

library(tidyverse)
dat <- data.frame(q=q, id_min = id_min, id_max=id_max)

p <- ggplot(dat, aes(x = q))

p + geom_line(aes(y = id_min), linetype=1) +
  geom_line(aes(y = id_max), linetype = 1) +
  ggthemes::theme_clean()

max(-1, min(p00 - (p01*p00/p11)-1,
            p11 - (p01*p10/p00)-1)
    )

#Q plausible values
q <- seq(p01,1-p10,.001)
# q = Pr(Z_star = 1) = True TGD Pop
q <- seq(0.006,0.03, .001)
id_min <- id_max <- rep(NA, length(q))

for (i in 1:length(q)) {
  id_min[i] <- max(-1*(p10 + p11)/(1-q[i]), 
                   -1*(p01/q[i])-(p10/(1-q[i])),
                   -1*(p00 + p01)/q[i])
  id_max[i] <- min((p00/(1-q[i]))-(p01/q[i]),
                   (p11/q[i])-(p10/(1-q[i])))
}

max(p00-((p01*p10)/p11)-1, p11-(p01*p10/p00)-1)

library(tidyverse)
dat <- data.frame(q=q, id_min = id_min, id_max=id_max)

#dat <- dat %>%
#  mutate(id_min=-1*id_min, id_max = -1*id_max)

p <- ggplot(dat, aes(x = q))

p + geom_line(aes(y = id_min), linetype=2) +
  geom_line(aes(y = id_max), linetype = 1) +
  geom_hline(yintercept = 0, color="red", linetype=2) +
  ggthemes::theme_clean() +
  ylim(-1,1)

#no measurement error
q <- p01+p11
id_min <- id_max <- rep(NA, length(q))

for (i in 1:length(q)) {
  id_min[i] <- max(-1*(p10 + p11)/(1-q[i]), 
                   -1*(p01/q[i])-(p10/(1-q[i])),
                   -1*(p00 + p01)/q[i])
  id_max[i] <- min((p00/(1-q[i]))-(p01/q[i]),
                   (p11/q[i])-(p10/(1-q[i])))
}

# lpsolve myself?
library(lpSolve)

const.mat <- matrix(c(1-.006, 0, .006, 0, 0,
                    0,0,1,1,0,
                    1,1,0,0,0,
                    0,1,1,0,1),
                    nrow = 4,byrow = TRUE)

const.dir <- c("=", "=", "=","=")

const.rhs <- c(p00, 1-(p01/.006), 1-(p10/(1-.006)), 1)

optimum <- lp(direction = "min", objective.in = c(1,0,0,1,-1),
              const.mat, const.dir,const.rhs)

optimum$solution

## as in the paper?
const.rhs <- c(p00, p01, p10,p11,0,0,1,
               0,0,0,0,0,0,0,0)
const.dir <- c("=","=","=","=", "=", "=","<",
               ">=",">=",">=",">=",">=",">=",">=",">=")
objective_in <- c(0, 0, -1, -1, 0, 0, 1, 1)

q <- seq(0.004,0.024, .0001)

id_min <- id_max <- rep(NA, length(q))

for (i in 1:length(q)){
  const.max <- matrix(c(1-q[i], 0, 0,0, q[i],0 ,0, 0,
                      0, 1-q[i], 0, 0, 0, q[i], 0, 0,
                      0, 0, 1-q[i], 0, 0, 0, q[i], 0,
                      0, 0, 0, 1-q[i], 0, 0, 0, q[i],
                      0, 1, 0, 0, 0, 0, 0, 0,
                      0, 0, 0, 0, 0, 0, 1, 0,
                      0, 1, 0, 1, 1, 0, 1, 0,
                      1, 0, 0, 0, 0, 0, 0, 0,
                      0, 1, 0, 0, 0, 0, 0, 0,
                      0, 0, 1, 0, 0, 0, 0, 0,
                      0, 0, 0, 1, 0, 0, 0, 0,
                      0, 0, 0, 0, 1, 0, 0, 0,
                      0, 0, 0, 0, 0, 1, 0, 0,
                      0, 0, 0, 0, 0, 0, 1, 0,
                      0, 0, 0, 0, 0, 0, 0, 1),
                    nrow = 15,
                    byrow=TRUE)

optimum <- lp(direction = "min", objective.in = objective_in,
              const.max, const.dir,const.rhs)
id_min[i] <- optimum$objval
optimum <- lp(direction = "max", objective.in = objective_in,
              const.max, const.dir,const.rhs)
id_max[i] <- optimum$objval
}

dat2 <- data.frame(q=q, id_min = id_min, id_max=id_max)
dat$id_max2 <- dat$id_max
dat$id_max2[dat$id_max<0] <- dat2$id_min[dat$id_max<0]

p <- ggplot(dat2, aes(x = q))

p +   geom_line(data=dat, aes(y=id_max2), linetype=1, color="blue") +
  geom_line(aes(y = id_min), linetype=2) +
  geom_line(aes(y = id_max), linetype = 1) +
  geom_hline(yintercept = 0, color="red", linetype=2) +
  annotate("point", x = 0.0148, y = .0877, color = "red") +
  ggthemes::theme_clean() +
  ylim(min(dat2$id_min),0.5) +
  xlim(0.002,0.025) +
  labs(x = expression(Treatment~Assignment~Probability~Q==Pr(Z[i]^{"*"}==1)),
       y = "Average Treatment Effect")

ggsave("auxplot1.tiff", device = "tiff", path = "C:/Users/aflor/Downloads/", width = 8, height = 6, unit = "in")


### Assumptions 1,3,4 ? 
const.rhs <- c(p00, p01, p10,p11,1,
               0,0,0,0)
const.dir <- c("=","=","=","=", "<",
               "=","=","=","=")
objective_in <- c(0, 0, 0, 0, 1, 1, 0, 0,
                  0, -1, 0, -1, 0, -1, 0, -1)

q <- seq(0.004,0.024, .0001)

id_min <- id_max <- rep(NA, length(q))

for (i in 1:length(q)){
  const.max <- matrix(c(0, 0, q[i], q[i], 0, 0, 0, 0, 1-q[i], 0, 1-q[i], 0, 0, 0, 0, 0,
                        0, 0, 0, 0, 0, 0, q[i], q[i], 0, 0, 0, 0, 1-q[i], 0, 1-q[i], 0,
                        q[i], q[i], 0, 0, 0, 0, 0, 0, 0, 1-q[i], 0, 1-q[i], 0, 0, 0, 0,
                        0, 0, 0, 0, q[i], q[i], 0, 0, 0, 0, 0, 0, 0, 1-q[i], 0, 1-q[i],
                        1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1,
                        1, rep(0,15),
                        0, 1, rep(0,14),
                        rep(0,12), 1, 0, 0, 0,
                        rep(0,14), 1, 0),
                      nrow = 9,
                      byrow=TRUE)
  
  optimum <- lp(direction = "min", objective.in = objective_in,
                const.max, const.dir,const.rhs)
  id_min[i] <- optimum$objval
  optimum <- lp(direction = "max", objective.in = objective_in,
                const.max, const.dir,const.rhs)
  id_max[i] <- optimum$objval
}

dat3 <- data.frame(q=q, id_min = id_min, id_max=id_max)

p <- ggplot(dat3, aes(x = q))

p + geom_line(aes(y = id_min), linetype=2) +
  geom_line(aes(y = id_max), linetype = 1) +
  geom_hline(yintercept = 0, color="red", linetype=2) +
  annotate("point", x = 0.0148, y = .0877, color = "red") +
  ggthemes::theme_clean() +
  ylim(min(dat2$id_min),0.5) +
  xlim(0.002,0.025) +
  labs(x = expression(Treatment~Assignment~Probability~Q==Pr(Z[i]^{"*"}==1)),
       y = "Average Treatment Effect")

### sensitivity approach

const.dir <- c("=","=","=","=", "=", "=","<",
               ">=",">=",">=",">=",">=",">=",">=",">=", "<=")
objective_in <- c(0, 0, -1, -1, 0, 0, 1, 1)

q <- 0.0148
rho <- seq(0,1,.001)

id_min <- id_max <- rep(NA, length(rho))

for (j in 1:length(rho)) {
  const.rhs <- c(p00, p01, p10,p11,0,0,1,
                 0,0,0,0,0,0,0,0, rho[j])
  for (i in 1:length(q)){
  const.max <- matrix(c(1-q[i], 0, 0,0, q[i],0 ,0, 0,
                        0, 1-q[i], 0, 0, 0, q[i], 0, 0,
                        0, 0, 1-q[i], 0, 0, 0, q[i], 0,
                        0, 0, 0, 1-q[i], 0, 0, 0, q[i],
                        0, 1, 0, 0, 0, 0, 0, 0,
                        0, 0, 0, 0, 0, 0, 1, 0,
                        0, 1, 0, 1, 1, 0, 1, 0,
                        1, 0, 0, 0, 0, 0, 0, 0,
                        0, 1, 0, 0, 0, 0, 0, 0,
                        0, 0, 1, 0, 0, 0, 0, 0,
                        0, 0, 0, 1, 0, 0, 0, 0,
                        0, 0, 0, 0, 1, 0, 0, 0,
                        0, 0, 0, 0, 0, 1, 0, 0,
                        0, 0, 0, 0, 0, 0, 1, 0,
                        0, 0, 0, 0, 0, 0, 0, 1,
                        0, 1, 0, 1, 1, 0, 1, 0),
                      nrow = 16,
                      byrow=TRUE)
  
  optimum <- lp(direction = "min", objective.in = objective_in,
                const.max, const.dir,const.rhs)
  id_min[j] <- optimum$objval
  optimum <- lp(direction = "max", objective.in = objective_in,
                const.max, const.dir,const.rhs)
  id_max[j] <- optimum$objval
}
}

dat4 <- data.frame(rho=rho, id_min = id_min, id_max=id_max)

p <- ggplot(dat4, aes(x = rho))

p +  geom_line(aes(y = id_min), linetype=2) +
  geom_line(aes(y = id_max), linetype = 1) +
  geom_hline(yintercept = 0, color="red", linetype=2) +
  annotate("point", x = 0, y = .0877, color = "red") +
  ggthemes::theme_clean() +
  xlim(0,.01) +
  labs(x = expression(rho),
       y = "Average Treatment Effect")

## maybe no assumption 5?
const.dir <- c("=","=","=","=","<",
               ">=",">=",">=",">=",">=",">=",">=",">=", "<=")
objective_in <- c(0, 0, -1, -1, 0, 0, 1, 1)

q <- 0.0148
rho <- seq(0,1,.001)

id_min <- id_max <- rep(NA, length(rho))

for (j in 1:length(rho)) {
  const.rhs <- c(p00, p01, p10,p11,1,
                 0,0,0,0,0,0,0,0, rho[j])
  for (i in 1:length(q)){
    const.max <- matrix(c(1-q[i], 0, 0,0, q[i],0 ,0, 0,
                          0, 1-q[i], 0, 0, 0, q[i], 0, 0,
                          0, 0, 1-q[i], 0, 0, 0, q[i], 0,
                          0, 0, 0, 1-q[i], 0, 0, 0, q[i],
                          0, 1, 0, 1, 1, 0, 1, 0,
                          1, 0, 0, 0, 0, 0, 0, 0,
                          0, 1, 0, 0, 0, 0, 0, 0,
                          0, 0, 1, 0, 0, 0, 0, 0,
                          0, 0, 0, 1, 0, 0, 0, 0,
                          0, 0, 0, 0, 1, 0, 0, 0,
                          0, 0, 0, 0, 0, 1, 0, 0,
                          0, 0, 0, 0, 0, 0, 1, 0,
                          0, 0, 0, 0, 0, 0, 0, 1,
                          0, 1, 0, 1, 1, 0, 1, 0),
                        nrow = 14,
                        byrow=TRUE)
    
    optimum <- lp(direction = "min", objective.in = objective_in,
                  const.max, const.dir,const.rhs)
    id_min[j] <- optimum$objval
    optimum <- lp(direction = "max", objective.in = objective_in,
                  const.max, const.dir,const.rhs)
    id_max[j] <- optimum$objval
  }
}

dat5 <- data.frame(rho=rho, id_min = id_min, id_max=id_max)

p <- ggplot(dat5, aes(x = rho))

p +  geom_line(aes(y = id_min), linetype=2) +
  geom_line(aes(y = id_max), linetype = 1) +
  geom_hline(yintercept = 0, color="red", linetype=2) +
  annotate("point", x = 0, y = .0877, color = "red") +
  ggthemes::theme_clean() +
  xlim(0,1) +
  labs(x = expression(rho),
       y = "Average Treatment Effect")

ggsave("auxplot2.tiff", device = "tiff", path = "C:/Users/aflor/Downloads/", width = 8, height = 6, unit = "in")

## sensitivity approach from the PSRM article as well?
# need to articualte the mis-measurement mechanism

cces <- haven::read_dta("C:/Users/aflor/Documents/CCES_2018/pooled_cces.dta") %>%
  filter(is.na(new_trans)==FALSE & is.na(prob_any)==FALSE)

mod_1 <- glm(prob_any ~ new_trans + factor(year), data = cces,
             family = binomial(link="logit"))
summary(mod_1)

mod_2 <- glm(prob_any ~ new_trans + log_age + factor(gender) + factor(race_2) + factor(educ_2) + factor(year),
             data = cces, family = binomial(link="logit"))
summary(mod_2)

library(MASS)
library(matrixStats)

etavec <- seq(0.001, 0.99, length.out = 20)

results <- NULL
results2 <- NULL

for (i in 1:length(etavec)){
  res <- NULL
  res2 <- NULL
  count <- 1
  while(count <= 250){
    # create data given eta: No knowledge of Ci.
    cces$newtreat <- NA
    cces$newtreat[cces$new_trans==1] <- rbinom(length(cces$new_trans[cces$new_trans==1]), 1, 1-etavec[i])
    cces$newtreat[cces$new_trans==0] <- cces$new_trans[cces$new_trans==0]  
    
    m <- glm(prob_any ~ newtreat + factor(year), family = binomial(link="logit"), data = cces)
    m2 <- glm(prob_any ~ newtreat + age + factor(gender) + factor(race_2) + factor(educ_2) + factor(year),
                 data = cces, family = binomial(link="logit"))
    res <- rbind(res, c(coef(m)["newtreat"], sqrt(vcov(m)[2,2])))
    res2 <- rbind(res2, c(coef(m2)["newtreat"], sqrt(vcov(m2)[2,2])))
    count <- count + 1
  }
  results <- rbind(results, c(etavec[i], colMedians(res)))
  results2 <- rbind(results2, c(etavec[i], colMedians(res2)))
}

colnames(results) <- colnames(results2) <- c("eta","est","se")
results <- as.data.frame(results)
results2 <- as.data.frame(results2)

results$lb <- results$est + qnorm(.025)*results$se
results$ub <- results$est + qnorm(.975)*results$se

results2$lb <- results2$est + qnorm(.025)*results2$se
results2$ub <- results2$est + qnorm(.975)*results2$se

p <- ggplot(results, aes(x = eta, y = est, ymin = lb, ymax = ub))

p1 <- p + geom_ribbon(alpha = 0.6, fill = "black") +
  geom_hline(yintercept=0, color = "red", linetype =2) +
  geom_line() + 
  ggthemes::theme_clean() +
  labs(x = expression(eta),
       y = "Coefficient") +
  ylim(-0.25, 3.25)

p <- ggplot(results2, aes(x = eta, y = est, ymin = lb, ymax = ub))

p2 <- p + geom_ribbon(alpha = 0.6, fill = "black") +
  geom_hline(yintercept=0, color = "red", linetype =2) +
  geom_line() + 
  ggthemes::theme_clean() +
  labs(x = expression(eta),
       y = "Coefficient") 

cowplot::plot_grid(p1, p2, ncol=2, labels = "auto")

ggsave("auxplot3.tiff", device = "tiff", path = "C:/Users/aflor/Downloads/", width = 8, height = 6, unit = "in")

### now to do the MLM maybe select limited

## Using the Stata results

dat_2 <- haven::read_dta("C:/Users/aflor/OneDrive - american.edu/Voting while Transgender/modpred_2.dta")

my_est <- NULL
my_est <- t(apply(dat_2, 2, quantile, c(0,.5,1)))
colnames(my_est) <- c("lb", "est","ub")
dat_2 <- as.data.frame(my_est)
dat_2$trans <- rep(c("Cisgender", "TGD"), 4)
dat_2$vid <- c(rep("Less Strict",4), rep("More Strict",4))
dat_2$id <- rep(c(rep("Easier to Update ID",2),rep("Difficult to Update ID",2)),2)

dat_2 <- dat_2 %>%
  mutate(id = factor(id, levels = c("Easier to Update ID", "Difficult to Update ID")))

p <- ggplot(dat_2, aes(x=vid, y = est, ymin = lb, ymax = ub, color = trans))

p1 <- p + geom_pointrange() +
  scale_color_brewer(type = "qual", palette=2) +
  ggthemes::theme_clean() +
  theme(legend.background = element_blank(), legend.position = "bottom") +
  labs(x = "Strictness of Voter ID Law",
       y = "Pr(Any Problem)",
       color = NULL) +
  facet_wrap(~id)

dat_3 <- haven::read_dta("C:/Users/aflor/OneDrive - american.edu/Voting while Transgender/modpred_5.dta")

my_est <- NULL
my_est <- t(apply(dat_3, 2, quantile, c(0,.5,1)))
colnames(my_est) <- c("lb", "est","ub")
dat_3 <- as.data.frame(my_est)
dat_3$trans <- rep(c("Cisgender", "TGD"), 4)
dat_3$vid <- c(rep("Less Strict",4), rep("More Strict",4))
dat_3$id <- rep(c(rep("Easier to Update ID",2),rep("Difficult to Update ID",2)),2)

dat_3 <- dat_3 %>%
  mutate(id = factor(id, levels = c("Easier to Update ID", "Difficult to Update ID")))

p <- ggplot(dat_3, aes(x=vid, y = est, ymin = lb, ymax = ub, color = trans))

p2 <- p + geom_pointrange() +
  scale_color_brewer(type = "qual", palette=2) +
  ggthemes::theme_clean() +
  theme(legend.background = element_blank(), legend.position = "bottom") +
  labs(x = "Strictness of Voter ID Law",
       y = "Pr(Any Problem)",
       color = NULL) +
  facet_wrap(~id)

dat_4 <- haven::read_dta("C:/Users/aflor/OneDrive - american.edu/Voting while Transgender/modpred_8.dta")

my_est <- NULL
my_est <- t(apply(dat_4, 2, quantile, c(0,.5,1)))
colnames(my_est) <- c("lb", "est","ub")
dat_4 <- as.data.frame(my_est)
dat_4$trans <- rep(c("Cisgender", "TGD"), 4)
dat_4$vid <- c(rep("Less Strict",4), rep("More Strict",4))
dat_4$id <- rep(c(rep("Easier to Update ID",2),rep("Difficult to Update ID",2)),2)

dat_4 <- dat_4 %>%
  mutate(id = factor(id, levels = c("Easier to Update ID", "Difficult to Update ID")))

p <- ggplot(dat_4, aes(x=vid, y = est, ymin = lb, ymax = ub, color = trans))

p3 <- p + geom_pointrange() +
  scale_color_brewer(type = "qual", palette=2) +
  ggthemes::theme_clean() +
  theme(legend.background = element_blank(), legend.position = "bottom") +
  labs(x = "Strictness of Voter ID Law",
       y = "Pr(Any Problem)",
       color = NULL) +
  facet_wrap(~id)

cowplot::plot_grid(p1,p2,p3, labels="auto", nrow=1, scale = 0.9)

ggsave("auxplot4.tiff", device = "tiff", path = "C:/Users/aflor/Downloads/", width = 12, height = 6, unit = "in")

